Mesoscopic Multi-Particle Collision Dynamics of Reaction-Diffusion Fronts 
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A mesoscopic multi-particle collision model for fluid dynamics is generalized to incorporate the 
chemical reactions among species that may diffuse at different rates. This generalization provides 
a means to simulate reaction-diffusion dynamics of complex reactive systems. The method is illus- 
trated by a study of cubic autocatalytic fronts. The mesoscopic scheme is able to reproduce the 
results of reaction-diffusion descriptions under conditions where the mean field equations are valid. 
The model is also able to incorporate the effects of molecular fluctuations on the reactive dynamics. 
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I. INTRODUCTION 

Mesoscopic models provide coarse-grained descriptions 
of the dynamics of systems that neglect certain details 
at microscopic scales while retaining essential dynamical 
features at mesoscopic and macroscopic scales. Conse- 
quently, a convenient way to study of the dynamics of 
complex systems over a large range of interesting space 
and time scales is through the use of such models. In 
physical and biological systems we often encounter situ- 
ations where mean field descriptions of reactions break 
down and molecular fluctuations play an important role 
in determining the character of the system's dynamics. 
Such effects are especially relevant for reactions taking 
place in nano-scalc domains or biochemical reactions at 
the cellular level. Fluctuations also play a role in far- 
from-equilibrium systems near bifurcation points or when 
the system behaves chaotically since the system is es- 
pecially susceptible to perturbations in such regimes. 1 
Mesoscopic models are able to capture the influence of 
such molecular fluctuations on the dynamics. Meso- 
scopic models are also useful for simulating the dynamics 
of macroscopic systems because they often provide sta- 
ble particle-based simulation schemes and can be imple- 
mented in complex geometries. 

In this article we consider a generalization of a meso- 
scopic multi-particle collision (MPC) (or stochastic rota- 
tion) mode&i4 to a pattern-forming chemically reacting 
system. We show how the multi-particle collision rule 
can be generalized to a multi-component system to yield 
different diffusion coefficients for the chemical species. 
Differences in diffusion coefficients can give rise to chem- 
ical instabilities which cannot occur if the diffusion coeffi- 
cients of all species are equal. Reactions are incorporated, 
also at a mesoscopic level, by combining a birth-death de- 
scription of reactive events with multi-particle collisions. 
The mesoscopic dynamics preserves all the basic conser- 
vation laws of the system and leads to the macroscopic 
evolution laws on long distance and time scales. 

To illustrate the scheme, the reactive MPC dynamics 
is used to investigate the evolution and structure of a cu- 



bic autocatalytic front. The cubic autoatalytic reaction 
is A + 2B — > 3B, where the autocatalyst B consumes the 
fuel A. If one considers a two-dimensional rectangular 
domain (or a thin rectangular slab in three dimensions) 
with B in left portion and A in the right portion, a re- 
action front will propagate from left to right. While the 
simulations presented in this paper are for cubic auto- 
catalytic fronts, the manner in which the diffusion pro- 
cess is modelled to yield different diffusion coefficients for 
different chemical species and the way reactions are in- 
corporated in the model presage extensions of the theory 
and applications to more complex far-from-equilibrium 
reactive systems. 

The paper is organized as follows: In Sec. [H] we 
sketch the basic elements of the multi-particle collision 
model and present its generalization to reactive systems 
where the chemical species can have different diffusion 
coefficients. Section IIIII describes the simulation of cu- 
bic autocatalytic fronts and compares the results of the 
mesoscopic simulations with the predictions of reaction- 
diffusion equations. The conclusions of the paper are 
given in Sec. IIVI 



II. MESOSCOPIC MODEL 

In multi-particle collision dynamics a system contain- 
ing N particles with continuous positions and velocities 
Vj evolves through a sequence of free streaming and colli- 
sion stepsi The collisions among the particles take place 
in the following way: the system is divided into cells and 
at time intervals r each cell labelled by £ is assigned at 
random a rotation operator cD^ from some suitable set of 
rotation operators. The center of mass velocity of 
the particles in cell £ is computed and the post-collision 
velocity of particle i in the cell is determined by rotat- 
ing its velocity, relative to the cell center of mass velocity, 
and adding the center of mass velocity to the result of this 
rotation: 



(1) 
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The velocity of every particle in cell £ is rotated by 
the same rotation operator but the rotation operator 
varies from cell to cell. The dynamics then consists 
free streaming interspersed by these multi-particle col- 
lision events. It has been shown that this dynamics con- 
serves mass, momentum and energy and thus leads to 
the full set of Navier-Stokes equations on long distance 
and time scales^^. The method has been applied to the 
study of a variety of systems^ including hydrodynamic 
flow3&, colloids^, polymers^, Brownian motion^ and sim- 
ple diffusion-influenced reaction dynamics^. 

We present a generalization of this model that allows 
the dynamics of reaction-diffusion systems to be investi- 
gated. This generalization entails several extensions of 
the MPC model. In particular, a multi-component ver- 
sion of the MPC modeliMi must be constructed that 
accounts for reactions among the chemical species and 
allows for the possibility that the diffusion coefficients of 
the species differ. 



Diffusion 
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Figure 1: Diffusion coefficient Da (7a) as a function of the 
density ua for various values of "/a- Squares, "/a ~ 1.00. 
Circles, "/a = 0.75. Triangles, 7a = 0.50. Diamonds ja = 
0.25. Solid lines plot the theoretical value (Eq. ©) of the 
diffusion coefficient. 



A multi-component MPC dynamics that provides a 
simple way to control the diffusion coefficients of different 
chemical species can be constructed as follows. Suppose 
we have s species labelled by an index a. Instead of ap- 
plying the MPC operator to all particles in a cell, we as- 
sume that multi-particle collision operators act to change 
the velocities of a fraction of the particles of species a in 
a cell for a — 1, . . . , s. More specifically, in each cell £ 
each particle of species a is chosen with probability 7 Q . 
If vf is the velocity of a chosen particle i of species a 
and V| is the center of mass velocity of all chosen par- 
ticles, the post-collision velocities of those particles that 
undergo collision are given by 



= Vf +<J> € (v?-V£) 



(2) 



The post-collision velocities of the particles that do not 
take part in the multi-particle collision do not change. 
The diffusion coefficients D a are functions of {7 Q '|o/ = 
1, . . . , s}, which can be tuned to change the values of the 
diffusion coefficients. 

In order to investigate the range over which the diffu- 
sion coefficients can vary, we consider the self diffusion 
coefficient of a single species A and change both the mean 
particle density ua and the fraction "/a of particles that 
participate in the multi-particle collisions. Figure^plots 
Da{ia), determined from the slope of the mean square 
displacement versus time, as a function density ua for 
different values of ja- /.From these results one sees that 
the self diffusion coefficient can be varied by about a fac- 
tor of five by changing the values of ja at a fixed density. 

The self diffusion coefficient for ja ^ 1 can be esti- 
mated in the Boltzmann approximation where correla- 
tions are neglected. The discrete-time Green-Kubo ex- 



pression for the diffusion coefficient is 4 *^ " 

Da^a) = M) + 52(v x (n)v x ) , (3) 



where, without loss of generality, we have set r = 1. Tak- 
ing into account the collision rule where, on average, a 
fraction ja of the particles undergo multi-particle colli- 
sions and fraction 1 — ja do not, we have 

Da{ia) = \(v 2 x ) + {1- 1 a)(v 2 x )+ia(v^v x ) + ... , (4) 

where t4 is the post-collision value of the velocity at 
time t = 1. Assuming that higher order collision terms 
can be expressed in terms of the first collision so the series 
is geometric, we obtain 

Da{ia) = -^(v 2 x ) + (v 2 x )(l + (l-j A )+ 7A r D + ...) 

(4) 
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(v x 1] v x ) 2(l-e- n *)+n A 
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3n A 



(5) 



(6) 



was computed in RefJS. The comparison in Fig. ^ shows 
that this analytical expression (solid lines) accurately de- 
scribes the simulation data. 
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Reaction 

The mesoscopic dynamics must also be generalized to 
allow for chemical reactions among the species. Our 
earlier study of diffusion-influenced reactions^ was re- 
stricted to a simple A + C ^ B + C reaction that occurs 
when the A or B particles collide with catalytic spheres 
C . Since we are now interested in reactions that occur 
among the mesoscopic particles, we instead use a birth- 
death stochastic law to describe the reactive event a 1 ! 12 . 

Here we restrict our considerations to the cubic auto- 
catalytic reaction A + 2B — > 3B. Independently in each 
cell we assume the reaction takes place with probabil- 
ity Pr = knAnB^B — 1), where n a is the number of 
molecules of species a in a cell. The reactive dynamics 
in a cell is described by the Markov chain^ 

P(n,i+l) = £V(n|n')P(n',f), (7) 

n' 

where n = (ua, ^b), P(n, t) is the probability that there 
are n particles in the cell at time t and the transition 
matrix W is given by 

W(n|n') = kn' A n' B (n' B - 1) V^-i^^+i 

+(1 - kn' A n' B {n' B - l))5 nA , n ' A 8 nB ,n' B • (8) 

The (discrete time) rate of change of the mean density of 
species a, n a (t) — ^ n n Q ,P(n, t), is 

n a (t + 1) - n a (t) =J2n a (w(n\n') - J n , n ,)p(n,i).(9) 

n,n' 

We assume that the MPC non-reactive collisions are suf- 
ficiently effective to maintain a local equilibrium Poisso- 
nian distribution in the cells so that 

P(n,t) & Pt(n A ;n A (t))P t (n B ;n B (t)) , (10) 

where the local Poisson distribution is Pg(n;n(t)) — 
e~ n ^n{t) n /n\. If we insert the local Poissonian approx- 
imation for P(n, t) in the right hand side of Eq. ® for 
a = A we obtain the discrete-time mean-field rate law, 

n A (t + 1) - n A {t) = -kn A (t)n 2 B (t) . (11) 

A similar equation can be derived for species B. Thus, 
the mass action rate law will describe the dynamics pro- 
vided diffusion is sufficiently rapid compared to reac- 
tion so that a local Poissonian distribution of particles 
is maintained during the evolution of the reactive sys- 
tem. In this limit the discrete-time rate law will closely 
approximate the continuous-time mass action rate law. 

After the reaction step, the particles free stream us- 
ing the post-collision values of the velocities, taking into 
account the boundary conditions of the system. Once 
all the particles have been moved, the time advances one 
unit and the multi-particle collision and reaction steps 
are applied again. This mesoscopic dynamics conserves 
the total mass, momentum and energy of the system. 



III. SIMULATION OF CHEMICAL FRONTS 

In this section we show that the mesoscopic MPC 
model can be used simulate the dynamics of cubic auto- 
catalytic fronts on macroscopic scales where comparisons 
with the predictions of reaction-diffusion equations can 
be made. Cubic autocatalytic fronts have been studied 
often in the context of a coupled pair of reaction-diffusion 
equations for the A and B species. 14 i 15 i 16 i 17 i 18 i 19 i 20 The 
particular focus of many of these studies was on the trans- 
verse front instability that occurs when the diffusion co- 
efficient of the fuel is sufficiently larger than that of the 
autocatalyst: at a critical value of the diffusion coefficient 
ratio an instability will develop and the planar front will 
become nonplanar and exhibit complex dynamics. 

Our investigations will be confined to a simpler case 
of a binary mixture undergoing the cubic autocatalytic 
reaction. For such a reacting mixture the relevant macro- 
scopic field variables are the total mass density p(r, t) = 
Pa + Pb, the local concentration c(r, t) = Pa/ Pi the cen- 
ter of mass velocity v(r, t) and the energy density e(r, t). 
For the isothermal cubic autocatalytic reaction with no 
net fluid flow so that v(r, t) = 0, and taking equal masses 
for the A and B species, the macroscopic equation for the 
number density of A is 21 , 

d 

—n A {r,t) = -kn A n 2 B +DV 2 n Al (12) 

where D is the mutual diffusion coefficient. The equation 
for ns(r,i) is not independent and follows from number 
conservation, h A + ub = no- 

Front profile 

The simulations of the reaction front using the MPC 
model were carried out in a rectangular prism with length 
I = 200 along x, width w = 200 along y and height h = 5 
units along z. The system was open along its length 
x, periodic boundary conditions were imposed in the 
y-direction, and bounce-back reflection boundary condi- 
tions were imposed on the top and bottom of the prism 
along z. In order to initiate a chemical front, A particles 
were distributed uniformly in the right side of the prism, 
(h A (x > 100) = 10, n B (x > 100) = 0), while B parti- 
cles were uniformly distributed in left side of the prism 
(n A (x < 100) = 0, n B (x < 100) = 10). The velocities 
were chosen from a Maxwell-Boltzmann distribution with 
reduced temperature fcgT = /3 _1 = 1/3. 

Starting from this initial condition a reaction front will 
develop as the autocatalyst B consumes the fuel A in 
the reaction. The front will move with velocity c and 
it is convenient to study the front dynamics in a frame 
moving with velocity c. Propagating fronts are depicted 
in Fig. [2 which shows the concentration field at a given 
time instant. The upper two panels plot the front for two 
values of the reaction rate constant fc and 7a = 7s = 1- 
We see that for k = 0.0005 the front profile is much 
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Figure 2: Concentration field at a given time instant for 
k = 0.0005 (top left panel) and k = 0.001 (top right panel). 
The system size is 200 x 200 x 5. Lower panels show the con- 
centration field for k = 0.0005 and 7 a = 0.25, jb = 1 (left), 
and k = 0.0005, 7a = 1, 7b = 0.25 (right). The structure of 
the reaction zone can be seen in these figures. 



thicker than that for k = 0.001. This dependence is 
in accord with predictions based on a reaction-diffusion 
description of the front as can be seen from the analysis 
given below. 

The structure of these planar fronts can investigated 
quantitatively by studying the front dynamics in a frame 
moving with the front velocity, £ = x — ct, and averaging 
the concentration profile over the width (along y) of the 
front, ua(0 — f dy riji(£,y). Figure|3]plots Jia(0 for the 
two values of k used in Fig. /.From this figure we see 
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Figure 3: Plot of the front profile wa(£) versus £ for different 
values of k: left panel, k — 0.0005; right panel, k — 0.001. The 
system size is 200 x 200 x 5. The continuous line represents 
the theoretical value obtained from Eq. I14H . 

that a well-defined propagating reaction front is obtained 
and the width of the front decreases as the reaction rate 
increases relative to the diffusion rate. 

The front shape and velocity can be determined from 
the reaction-diffusion equation. For a planar front prop- 
agating along the x-direction, in a frame moving with the 



front velocity, the reaction-diffusion equation l|12|) is, 

(13) 

The front profile can be found analytically from the so- 
lution of this equation and is given by^S 

n A (0 = n Q (l + e-^ D y\ (14) 

where the front speed is c = {Dkn^/2) 1 / 2 . The profile 
for species B can be found from the conservation con- 
dition fiA + Tlb = no. Figure 01 compares this analyt- 
ical prediction with the simulation results of the MPC 
reaction-diffusion dynamics. For 7^4 =75 = 1 the mu- 
tual diffusion coefficient D is given by Eq. JSJ. There is 
good agreement between the simulation and analytical 
values for small k where the conditions for the validity 
of the mean field approximation are satisfied. For larger 
k values, such as k — 0.001 in the right panel of the fig- 
ure we see that there are deviations from the mean field 
result. For this value of k, the reaction is fast and there 
is a breakdown of the local Poissonian equilibrium and a 
reaction-diffusion description is not applicable. A similar 
breakdown is observed for very small k, for example for 
k = 0.0002, due to the fact that very few reactive events 
occur in the reaction zone of the front and fluctuations 
are important. 

The front velocity was determined from the simulation 
data as a function of k. In Fig.0]we plot the front velocity 
c versus k and compare the simulation results with the 
prediction c = {Dkn^/2) 1 / 2 . The front velocity agrees 
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Figure 4: (filled circles) Front velocity c as a function of k. 
(solid line) Front velocity from Eq. 11411 . 

with the simulation results for k < 0.003, although the 
front profile deviates slightly from the predicted value for 
somewhat smaller values of k (k < 0.001). 
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More microscopic aspects of the front structure and 
dynamics that are captured by the MPC model are il- 
lustrated in the lower two panels of Fig. [21 These figures 
plot snapshots of the front for k = 0.0005, the same value 
of k as in the top left panel of the figure, but for two 
different pairs of j a values, (7,4 = 0.25,7b = 1.0) and 
(7,4 = 1.0,7s = 0.25). Comparison of the lower panels 
of the figures, and also with the upper left panel, shows 
that the structures of the interfacial zones are different. 
In the MPC dynamics employed here, the diffusion of the 
species depends on their density and 7 Q . Since the den- 
sity of the species changes significantly in the interfacial 
zone, it is likely that a concentration-dependent mutual 
diffusion coefficient is required to describe this structure. 

IV. CONCLUSION 

The generalizations of the multi-particle collision 
model described here, and its extensions, allow one to 
study a variety of phenomena at the mesoscopic level. In 
particular, the ability to simulate the dynamics of multi- 
component systems whose diffusion coefficients can be 
different means that diffusion-driven instabilities, such as 
the transverse cubic autocatalytic front instability con- 
sidered in this paper, can be investigated. Since the 
mesoscopic MPC model preserves the basic conservation 



laws in the system, to study such instabilities requires 
the presence of a third solvent species so that there are 
two independent diffusion coefficients in the system. The 
method could also be used to study reactive and non- 
reactive binary fluid flows which also show interesting 
instabilities where fluctuations play a role near the onset 
of instabilities. 

The cubic autocatalytic reaction is simply one exam- 
ple of a much broader class of reaction-diffusion systems 
that can be studied using reactive versions of the meso- 
scopic multi-particle collision dynamics. In particular, 
more general reaction-diffusion dynamics in specific ge- 
ometries relevant for the materials science and biological 
applications may be carried out. The presence of flows 
can also be treated easily in this context. 

While we have focused primarily on parameter do- 
mains where mean field approximations are largely ap- 
plicable, one of the most interesting applications of the 
methodology introduced in this paper is to systems on 
mesoscales where particle numbers are small so that fluc- 
tuations play a crucial role in the dynamics and system 
geometry is important. 
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